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Abstract 

Using basic thermodynamic principles we derive a Cahn-Hilliard-Darcy model for 
tumour growth including nutrient diffusion, chemotaxis, active transport, adhesion, 
apoptosis and proliferation. The model generalises earlier models and in particular 
includes active transport mechanisms which ensure thermodynamic consistency. We 
perform a formally matched asymptotic expansion and develop several sharp interface 
models. Some of them are classical and some are new which for example include a 
jump in the nutrient density at the interface. A linear stability analysis for a growing 
nucleus is performed and in particular the role of the new active transport term is 
analysed. Numerical computations are performed to study the influence of the active 
transport term for specific growth scenarios. 

Key words. Tumour growth, diffuse interface model, Cahn-Hilliard equation, chemo¬ 
taxis, Darcy’s flow, matched asymptotic expansions, stability analysis, hnite element com¬ 
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1 Introduction 

In the last decades the understanding of tumour related illnesses has undergone a swift 
development. Nowadays tumour therapy can be adapted to the genetic fingerprint of the 
tumour, resulting in a “targeted therapy” that has dramatically improved the prognosis of 
many illnesses. While some important mutations in tumour genomes have been identified 
and exploited by modern tumour drugs, basic growth behaviours of tumours are still far 
from being understood, e.g. angiogenesis and the formation of metastases. The complexity 
of oncology has also attracted increasing interest of mathematicians, who are trying to 
find the appropriate equations to provide additional insights in certain aspects of tumour 
growth, see for example [6] and m- In this paper we want to introduce a new diffuse 
interface model for tumour growth, and compare the resulting system of partial differential 
equations to some other recent contributions [laiiaEaiMiEHiEniEaETiiMiiiaiiii]. 
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In order to obtain a tractable system of partial differential equations, we will in this 
paper neglect some effects which could be addressed in further research and which then 
would lead to more complete theories. From a medical point of view we will hence make 
the following assumptions as foundations for our further considerations; 

1. Tumour cells only die by apoptosis. Hence we neglect the possibility of tumour 
necrosis, where we would have to take account of the negative effects of chemical 
species from the former intracellular space on the surrounding tumour cells. 

2. The tissue around the tumour does not react to the tumour cells in any active way. 
In particular, we neglect any response of the immune system to the tumour tissue. 

3. Larger tumour entities are actually enforcing blood vessel growth towards themselves 
by secreting vessel growth factors. This is a phenomenon that could be addressed in 
future in a generalised model. 

4. We postulate the existence of an unspecified chemical species acting as a nutrient 
for the tumour cells. This nutrient is not consumed by the healthy tissue. We will 
also introduce terms which will reflect chemotaxis, which is the active movement 
of the tumour colony towards nutrient sources. Additionally, the introduction of 
chemotaxis will also lead to the opposite process, meaning that the nutrient is moving 
towards the nearby tumour cells. As we will point out later, this could be seen as a 
correlate of a nutrient uptake mechanism. 

Here we state a slightly simplified version of the general system, which will be derived 
in Sectionfrom thermodynamic principles. We will derive and analyse a two-component 
mixture model of tumour and healthy cells, whose behaviour is governed by the system 


divu = ar, (1.1a) 

u =-Ar(Vp-/uV(^-X<^cTVv5), (1-lb) 

dtif + div {v(p) = V • + psT, (1-lc) 

p= - l3sAip-x>f<T, (1-ld) 

dtcr + div (av) = div {n{Lp){xa ^cf - Xip^V>)) “ Cah{(p), (1-le) 

T = {ra-A)h{ip). (l.lf) 


Here, v denotes the volume-averaged velocity of the mixture, p denotes the pressure, a 
denotes the concentration of an unspecified chemical species that serves as a nutrient 
for the tumour, ip e [-1,1] denotes the difference in volume fractions, with {ip = 1} 
representing unmixed tumour tissue, and {p = -1} representing the surrounding healthy 
tissue, and p denotes the chemical potential for p. The particular simple form of (1.1a) 
is different to earlier modelling attempts and is based on the fact that we use volume- 
averaged velocities. 

The positive constants K, /3, 'P, A, and C denote the permeability, surface tension, 
proliferation rate, apoptosis rate, and consumption rate, respectively. The constants pg 
and a are related to the densities of the two components (see (2.32) below), in particular, 
for the case of matched densities we have a = 0. Meanwhile m{p) and n(p) are non¬ 
negative mobilities for p and a, respectively, and T(-) is a potential with two equal minima 
at ±1. In addition, we choose h as an interpolation function with h{-l) = 0 and h{\) = 1. 
The simplest choice is given as h{p) = ^{1 + p). 

We denote Xo- > 0 as the diffusivity of the nutrient, and Xv^ ^ 0 can be seen as a 
parameter for transport mechanisms such as chemotaxis and active uptake (see below for 
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more details). Finally, the parameter e is related to the thickness of the interfacial layers 


present in phase field systems. The system (1.1) is a Cahn-Hilliard-Darcy system coupled 


to a convection-diffusion-reaction equation for the nutrient. 

(1.1a) and (1.1b) model the mass balance using a Darcy-type system, and in the 
situation of unmatched densities (a t 0), the gain and loss of volume resulting from 
the mass transition F leads to sources and sinks in the mass balance. In (1.1c) and 


(l.ld), ip is governed by a Cahn-Hilliard type equation with additional source terms. 


The mass transition from the the healthy cells to the tumour component and vice versa is 


described in (l.lf), where tumour growth/proliferation is represented by the term Vah{ip), 


and the process of apoptosis is modelled by the term Ah((p). In (l.le), the nutrient 


is subjected to an equation of convection-reaction-diffusion type, and the term Cah{(p) 
represents consumption of the nutrient only in the presence of the tumour cells. As in 
[To] , we could also consider the situation where the tumour possesses its own vasculature 
and the nutrient may be supplied to the tumour via a capillary network at a rate 
where as is the constant nutrient concentration in the vasculature and B is the blood- 
tissue transfer rate which might depend on (p and x. This leads to the following nutrient 


balance equation instead of (l.le) 


dt(T + div(cru) = div(n((/?)(xcrVcr - -Cah{(p) +B{aB -cr). 


Under appropriate boundary conditions the system (1.1) allows for an energy inequality 


see (2.27) below) and we believe that this inequality will allow the well-posedness of the 


above system to be rigorously shown. 

We now motivate the particular choices for the modelling of proliferation, apoptosis, 


chemotaxis, and mass transition in (1.1) 


In (l.lf), we obtain that F = Va-A holds in the tumour region {</? = 1}. The implicit 


assumption that the tumour growth is proportional to the nutrient supply can be 
justihed by the fact that malign tumours have the common genetic feature that 
certain growth inhibiting proteins have been switched off by mutations. Hence, we 
can assume that while in healthy cells the mitotic cycle is rather strictly inhibited, 
tumour cells often show unregulated growth behaviour which is only limited by the 
supply of nutrients. 

Moreover, implicit in the choice of zero mass transition F = 0 in the healthy region 
{(p = -1} is the assumption that the tumour proliferation rate is more significant 
than that of the healthy tissue. 


In (I.lc) and (l.le), the fluxes for ip and a are given by 


■■= -m{ip)Vp = -m{ip)V (f - /3eA(y9 - 
Qa ■■= -n{ip)S/{XaCr - Xv<^), 


respectively. It has been pointed out by Roussos, Condeelis and Patsialou in uni 
that the undersupply of nutrient induces chemotaxis in certain tumour entities. This 
is reflected in the term m{ip)V{xifi^) of q>fi which drives the cells towards regions of 
high nutrient. 

On the other hand, we note that the term n{p)\I{xipp) in drives the nutrient 
to regions of high (p, i.e., to the tumour cells, which indicates that the nutrient is 
actively moving towards the tumour cells. This may seem to be counter-intuitive at 
first glance. However, this term will only contribute to the equation signihcantly in 
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the vicinity of the interface between the tumour and healthy cells. This allows the 
interpretation that the term n{ip)V(Xip^) reflects active transport mechanisms which 
move the nutrient into the tumour colony. Here we use the term “active transport” 
in the biological sense in order to indicate that some kind of mechanism is needed to 
maintain the transport (in contrast to passive transporters which are driven only by 
the concentration gradient of the substance). The additional mechanism allows cells 
to establish persisting concentration differences between different compartments. In 
particular we can expect that tumours, which have these active transporters on their 
cell membrane, are not dependent on diffusion but can establish high concentration 
of the vital nutrient even against the nutrient concentration gradient. 


Here we briefly give an example, where mechanisms like this have already been ob¬ 
served: Malign tumour cells often have a significantly increased need for glucose, a fact 
that is sometimes referred to as the Warburg effect. As a consequence of several mutations 
in the tumour genome, these cells can adapt to their high rate of glucose consumption in 
several ways. Apart from angiogenesis, which leads to a well perfused tumour environ¬ 
ment providing large amount of glucose, the tumour cells can also express (i.e., build) 
more glucose transporters, which provide an improved glucose transport through the cell 
membrane. Recently, both passive glucose transporters, so-called GLUT proteins, and 
active glucose transporters called SGLTs, have been observed on the cell membrane of 
several tumour entities. For a more detailed description regarding the GLUT transporters 
we refer to m, whereas SGLT expression of tumours has been described by |31j and |41j . 
In the system we will derive in this paper, the existence of passive nutrient transporters 
like the GLUTs is implicitly assumed by including nutrient diffusion. Apart from that, 
it will become more obvious in the corresponding sharp interface system that the term 
n{(p)V represents an active nutrient transport towards the tumour. 

We note that in (1.1), the mechanism of chemotaxis and active transport are connected 


via the parameter Xp- In principle, it is possible to decouple the two mechanisms. In order 
to do so, we introduce the following choice for the mobility n(ip) and diffusion coefficient 
Xa (see also Section 3.3.3 below): For A > 0 and a non-negative mobility V{ip), we set 


n{ip) = XV{if)Xp , Xa = X^Xp- 
Then, the corresponding fluxes for ip and a are now given as 

Qp = -m{p>)V (f ^''(y?) - /3eA(/? - Xpcr) , 
Qa = -T>{(p)V{a - Xif). 


( 1 . 2 ) 

(1.3a) 

(1.3b) 


For this choice, we can switch off the effects of active transport by sending A ->• 0, while 
preserving the effects of chemotaxis. 

We now compare the new model (1.1) and some of the previous diffuse interface models 
in the literature: 


• An important difference to other models is that we take a volume-averaged velocity 
which leads to the simple form divu = ccF, as a consequence of mass balance. In 
other models, this equation has to be replaced by a more complicated transport 
equation. 


Another significant difference is the presence of the term - div (n{(p)xpX/'p) in (l-le), 
which only appears in cases where chemotaxis and active transport are taken into 
account. As we have pointed out before, it represents active nutrient transport 
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towards the tumour. The corresponding nutrient equations in [HI O [23 [2S1 Eni 
(381 m] do not include an equivalent term. However, we point out that this active 
transport mechanism is present in the nutrient equation of [ 22 ]) ''^tio however used 
different source terms and no Darcy-flow contributions. 


Our choice of the mass transition term T in (l.lf) can also be found in 
Alternatively, one may consider equations of the form 


dt(p= div(m(v7)V/u)+ P(vp)(cr-xy?-//), 

dta = div (n((^)( Vcr - xVy?)) - P{^){cr - - //), 


where the chemical potential enters as a source term for the equations of y? and 
a. Here x ^ 0 is a constant, and P{-) denotes a non-negative proliferation function. 
This type of mass transition term appears in [29] and in naEaEO] with X = 0 - 


The presence of chemotaxis, represented by the term -XifO' in (l.ld) can also be found 
in the models of |14[ 128112^ I38j . while the corresponding Cahn-Hilliard systems in 
USES] [30113211371 Ej do not include an equivalent term. 


In |32l 1371 E] , the nutrient does not enter the Darcy law for v like in (1.1b). 


In the diffuse interface model (1.1), the parameter e is related to the thickness of the 


interfacial layer, which separates the tumour cell regions {(^ = 1 } and the healthy cell 
regions {(p = -1}. Hence, it is natural to ask if a sharp interface description of the problem 
will emerge in the limit e ^ 0. This means in the limit the interface between the tumour 
cells and the healthy cells is represented by a hypersurface of zero thickness. 

For convenience, suppose we take the mobilities m{(p) = mo, n{(p) = no to be constant. 
A formally matched asymptotic analysis will yield the following sharp interface limit from 
(1.1) (see Section]^ for more details): Let fix and Qh denote the tumour cell region and 
the healthy cell region, respectively, which are separated by an interface S. Then it holds 
that 


V = -KVp in Qx 


divu 


-moA/i 


dtcr + div (av) - reoXo-^'^ 
[v]h-u = 0, 

[/^]h = 0) 

2(-V + V ■ v) 

2^{-V + v-u) 
X<T 


{a{Va-A) in Hr, 
jo in Hh, 

Upg-a){Vao-A) in Hr, 

|o in Q.H, 

( -Ca in Hr, 

0 in Hh, 

[cr]^ = 2^, [p]^ = /37 k on S, 

rp 

mo [^Po]h' ^ on S, 

rp 

no ■ 1 / on S. 


(1.4a) 

(1.4b) 


(1.4c) 


(1.4d) 


(1.4e) 

(1.4f) 

(l-4g) 

(1.4h) 


Here, 7 is a constant related to the potential 'L (see (3.17) below), V denotes the normal 

rp I I 

velocity of S, n is the mean curvature of S, [f]jj denotes the jump of / from Hr to Qh 
across S (see (3. 16])), and u is the outward unit normal of S, pointing towards Hr. 
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In comparison to the formal sharp interface limits of mEnmi!, the most significant 
difference is the jump condition (1.4e)2. Let us remark on its physical meaning. Let (Tt 
and ajj denote the limiting values of the nutrient on the interface S from tumour cell 


regions and from the healthy cell regions, respectively. Then, (1.4e)2 implies that 


o 

(JT = O' H + 4-. 


Thus, if Xip is positive, then (1.4e)2 tells us that the tumour cells will experience a higher 
level of nutrient concentration than the healthy cells on the interface, which reflects the 


effect of the active transport mechanism in (l.le), attracting nutrients from the healthy 
cell regions into the tumour. 


If we consider the fluxes (1.3) in (O), then one obtains the sharp interface model ( |1.4[ ) 
with the following modification (see Section 3.3.3 for more details): Instead of (1.4e)2, we 
now have 


\o^H - 2 ^- 


In particular, the parameter A only enters explicitly in the jump condition for a, which 
relates to the above discussion regarding the physical interpretation of (1.4e)2. Note 


that A ^ 0 is a consequence of the active transport term we have discussed in the phase 
field model. In the matched asymptotics expansion, this term directly leads to a jump 
of the nutrient concentration at the tumour interface. Therefore, we obtain exactly the 
situation one would expect from active transport mechanisms: Close to the tumour surface, 
we observe a higher nutrient concentration inside the tumour than on the outside of the 
tumour, a situation that is only possible due to the transporter molecules. Hence it should 
be considered if A could be referred to as a density parameter for the active transport 
proteins. 

The plan of this paper is as follows: In Section we derive the new phase field model 
from thermodynamic principles and compare with previous phase field models of tumour 
growth in the literature. In Section]^ we perform a formal asymptotic analysis to derive 
certain sharp interface models of tumour growth. In Section we investigate the stability 
of radial solutions to a particular sharp interface model via a linear stability analysis, 
and highlight the effect of the active transport parameter on the stability. In Section 
we present quantitative simulations for radially symmetric solutions and qualitative 
simulations for more general scenarios. 


2 Model Derivation 


Let us consider a two component mixture consisting of tumour and healthy cells in a 
bounded domain d= 1,2,3. 

We denote the first component as the component of healthy tissues, and the second 
component as the tumour tissues. Let pi, i = 1,2, denote the actual mass of the component 
matter per volume in the mixture, and let pi, i = 1,2, be the mass density of a pure 
component i. Then, p '■= pi + p 2 denotes the mixture density (which is not necessarily 
constant), and we define the volume fraction of component i as 


Ui = 


fH 

Pi 


( 2 . 1 ) 


We expect that physically, pi e [0,/9i] and thus Ui e [0,1]. In addition to the considerations 
stated in Section we make the following modelling assumptions: 
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• There is no external volume compartment besides the two components, i.e., 

ui+U2 = 1. (2-2) 

• We allow for mass exchange between the two components. Growth of the tumour 
is represented by mass transfer from component 1 (healthy tissues) to component 2 
(tumour tissues), while tumour cells are converted back into the surrounding healthy 
tissues when they die. 

• We choose the mixture velocity to be the volume-averaged velocity: 


V ■■= uivi + U 2 V 2 , (2-3) 

where Vi is the individual velocity of component i. 

• We model a general chemical species which is treated as a nutrient for the tumour 
tissues. Its concentration is denoted by a and it is transported by the volume- 
averaged mixture velocity and a flux J^-. 


2.1 Balance laws 

The balance law for mass of each component reads as 


dtpi + div(piui) = Ti, 

dtP2 + div(p2^2) = r2. 


Observe that by (2.1), we can write (2.4) in the following way: For i = 1,2, 

Fj 

dtUi + dw{uiVi) = 

Pi 


We see that by (2.2), (2.3), and (2.5), 

divu = div(uiui) + div {U 2 V 2 ) = ^ f ^ - dtUi] = 

i=i V Pi > 

We introduce the fluxes: 


r2 Fi ^ 
— + — F„. 

P2 Pi 


Ji Piip^i 1 ^ Jl 'i' J 2 '! J — _ + _ T2. 

Pi P 2 


Then, we see that 


(2.4a) 

(2.4b) 


(2.5) 


( 2 . 6 ) 


(2.7) 


J + pv = Ji + J 2 + pv = piVi + P2V2, 


and so, upon adding the equations in (2.4) we obtain the equation for the mixture density: 

dtp+ div{pivi + P 2 V 2 ) = dtp+ div{pv + J') = Ti - 1 -F 2 . (2.8) 

We now want to derive an equation for the phase field variable p. Recalling pi = piUi, we 
obtain from (2.5) that 


1 

dtUi + — div Jj -H div (uiv) = 

pi pi 


(2.9) 
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We define the order parameter ip as the difference in volume fractions: 


p-.= U2-U1, 


( 2 . 10 ) 


then, subtracting the equation for ui from the equation for M 2 , and using (2.7), we obtain 
the equation for p-. 


r2 Ti 

dtp + div((/ 9 u) + div J =3 -— =: 

P2 Pi 


( 2 . 11 ) 


We point out that from the constraint (2.2), we obtain 

1 + p 1-p 

"“2 = —ui = ——. 


Thus, the region of the tumour tissues is represented by {x e = 1} and the region of 
healthy tissues is represented by {x e : (^ = -1}. In particular, the mixture density p can 
be expressed as a linear function of p\ 


P = Pi + P2 = Pi- 


1 


P _ 1 + p 
— +P2^— 


Pi + P2 P2- Pi 
—;;— + —;;—■ 


( 2 . 12 ) 


For the nutrient, we postulate the following balance law: 


+ div((Tu) + div Jo-=-5, (2.13) 

where S denotes a source/sink term for the nutrient. In addition, av models the transport 
by the volume-averaged velocity and models other transport mechanisms like diffusion 
and chemotaxis. 


2.2 Energy inequality 

We postulate a general energy density of the form: 

e{p,Vp,a) = f{p,Vp) + N{p,a). (2.14) 


Here, we neglected inertia effects, and so the kinetic energy does not appear in e. Instead 
we refer the reader to [l3] for the derivation of a model that includes inertia effects, leading 


to a Navier-Stokes-Cahn-Hilliard version of (1.1). The first term / in (2.14) accounts for 
interfacial energy and unmixing tendencies, while the second term N describes the chemical 
energy of the nutrient and energy contributions resulting from the interactions between 
the tumour tissues and the nutrient. The latter will, for example, lead to chemotatic 
effects which are of particular interest as they result in the tumour tissue growing towards 
regions with high nutrient concentration. 

In the following, we will consider / to be of Ginzburg-Landau type: For constants 
A,B > 0, we choose 


f{p,\/p)--=A^{p) + ^\Vpf, (2.15) 

where 4'(s) is a potential with equal minima at s = ± 1 . 

We will now derive the diffuse interface model based on a dissipation inequality for 
balance laws with source terms which has been used similarly by Gurtin |25l [26] and Podio- 
Guidugli [39] to derive phase field and Cahn-Hilliard type equations. These authors used 
the second law of thermodynamics which in an isothermal situation is formulated as a 
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free energy inequality. We also refer to Chapter 62 of Curtin, Fried, and Anand \27\ for 
a detailed discussion of situations with source terms. The second law of thermodynamics 
in the isothermal situation requires that for all volumes V{t) c n which are transported 
with the fluid velocity the following inequality has to hold (see [251 Eg EH [39]): 

f edx<- r Je-i/d7i^~^+ f + C5(-5)dx, 

dt Jv(t) Jav(t) Jv{t) 

where v is the outer unit normal to dV{t) and Je is an energy flux yet to be specified. 
Following 1211, we postulate that the source terms F.^, F^, and the nutrient supply (-5) 
carry with them a supply of energy described by 

J^^^^CvT^ + CipT^ + cs{-S)d:^, (2-16) 

for some c.v,c^p and cs yet to be determined. 

Using the transport theorem and the divergence theorem, we obtain the following local 
form 


dte + div (eu) + div Je - c„F„ - CtpT^p + csS < 0. 


(2.17) 


We now use the Lagrange multiplier method of Liu and Muller, see for example Section 2.2 
of |2] and Chapter 7 of |36| . Let A„, Ao- and A<^ denote Lagrange multipliers for the diver¬ 
gence equation (2.6), the nutrient equation (2.13) and the order parameter equation (2.11). 
We require that the following inequality holds for arbitrary ((/J,c^,^),F^,,F(^,^S,0*y?,9*cT): 

-V ■■= dte + div (eu) + div Je - CvT^ - + csS 

- \v{ divu - F.u) 

- Xcr{d*a + fidivu -i- div J^ + S) 

- Xtp{d*ip + {pd\\v + div J - F,^) < 0, (2.18) 

where we used the notation 

dlip.= dt(p + Vip-v, 

as the material derivative of (p with respect to v. Using the identity 

Vip ■ d*{Vp) = div {d*pVp) - d*pdiv (Vp) - {Vp 0 Vp) ■ Vv, 
we compute that 

-V = div (Je - XipJ - Ao-Jcr + Bd*pVp) 

/ ON \ 

+ |^A'I''(y?) + — - BAp - Xipj d*p -Vv ■■ B {SJp 0 Vp) 

IdN \ 

+ I ~ Aq- j C + S Cg ~ A(j ) + V X^p ' J + V Act ' Jcr 

+ (e Xpp X(j(7 A^) div u + F.y (A.|; ) + F,^(A(^ Cp ). (2.19) 

We use the following notation: 

dN dN 

{p) + BAp. 
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Applying the product rule to the divergence term in (2.19), we then obtain 


-V = div (Je - XtpJ - Act Jo- + BdlipVif + (e - X^pip - X^a - A„)u) 

+ (^ ~ Ac^) (p S(^cg — Act) + r.|;(A.|; ~ Ctj) + ~ (-A^iO- “ Ao-) o' 

- Vu : B {\/ip 0 \/(p) - u ■ V(e - Xpp> - XaO - X^) + VX^ ■ J + VAo- ■ Ja- (2.20) 
Employing the following identities 


Vr’:(V</J0V<y9)= div ((V(^ 0 \/p>)v) - v ■ div (Vy? 0 Vy?), 

-V (IVvjI^) = div(V(^ 0 Vip) - A(pV(p, 

dt(pV(p = d*(pV<p - (Vcp ■ v)V(p = d*p>Vp> - (Vy? 0 V^p)v, 


in (2.20), we arrive at 


-V = div (Je - XpJ - Act Jo- + BdtipVip + (e - Xp(p - X^o - Xv)v) 

+ (^ ~ Ai^) (p + S(^cg — Act) + r.CT(A.CT ~ C'u') + r(^( Ai^ ~ (-A^,o- “ Ao) o 

- V • (v(e - Xpp> - ActCt - A^ - f |V<y9|^) - BAp>\lLp) + VA<^ • J -+ VAo • Jo- (2.21) 


2.3 Constitutive assumptions and the general model 


We are now seeking for a model fulfilling the second law of thermodynamics in the version 
of a dissipation inequality stated in Section 2.2 We do not aim for the most general model 
but will state certain constitutive assumptions which take the most relevant effects into 
account. We hence make the following constitutive assumptions: 


Je = XpJ + ActJct - BdtipVp) - (e - X^pip - X^a - Xv)v, (2.22a) 

= Act — A^,ct? Op = Xp = //, Cu = A-ct, (2.22b) 

Jo- = -n{(p)VN^o-, J = -m{ip)Vp, (2.22c) 


where n(y?) and m{ip) are non-negative mobilities. We introduce a pressure-like function 
p and choose 


Xv=p- A'^{p) - ^ + e- pip- N^„o, (2.23) 

and for a positive constant K, 

V = K {y{e- pp- - A„ - ^ - BApVp) 

= K (V(-p + A'^{p)) - BApS/p) 

= -K{\/p-{p-N^p)Vp). (2.24) 


Eq. (2.22a) makes a constitutive assumption for the energy flux Jg which guarantees that 
the divergence term in (2.21) vanishes. It contains classical terms like pJ and A^^o-Jo- which 
describe energy flux due to mass diffusion and the non-classical term BdtpVp which is 
due to moving phase boundaries, see also mm where this term is discussed. The last term 
in (2.22a) will result in energy changes due to work by macroscopic stress, compare [2]. 
Meanwhile, (2.22b), (2.22c), (2.23) and (2.24) are considered in order for the right hand 
side of (2.21) to be non-positive for arbitrary values of {p,a,v,rv,Tp,S,d*p,d*a). We 
mention that (2.24) is a Darcy law with force {p - N^p)'\/p. 
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Thus, the model equations for tumour growth are 


divu = r^, (2.25a) 

V =-K{Vp- (2.25b) 
dtp + div (pv) = div {m(p)V p) + T^,, (2.25c) 

p = A'I>'(p) - BAp + N^ip, (2.25d) 

dt(7 + div (cru) = div (n((^)VA^,o-) -S, (2.25e) 

where 

+ P2^r2, ^ip = p2^^2- 
Supplemented with the boundary conditions 

\/p ■ v = Vp ■ 1 / = 0 on (2.26) 


then the above model satishes the following energy equality: 
d 


dt Jn 


f 

Jn 


A^{p) + ^\Vpf + N{p,a) 


dx 


m{p)\Vp\^ + n{p) I ViV,, 


|2 1 I |2 1 

I H-|u| dx + 

K 


- Xv'^v - p'^^dx 


f V ■ u{N{p,a) + p) - n{p)N^^SjN^a ■ v ^ = 0. 
J d^l ’ ’ 


(2.27) 


This follows from integrating (2.21) over Vl and using the dehnition of -V from (2.18), the 
constitutive assumptions (2.22), (2.23), and (2.24), and applying the divergence theorem. 
Here, we have not prescribed boundary conditions for N^a- and v. We will look at suitable 
boundary conditions for them later. 

We point out that using (2.8), (2.12), (2.25a), (2.25c), and the definition of T,^ and 
r„, we obtain 

Ti + r 2 = dtp + div {pv) + div JT” 


———{dip + p div u) + ——— div v + div J' 
^2 “ Pi / J- , / / X , --It 


2 

div i^J 


{d\\ {m{p)Vp) + P2^P2 - /o/bi) + 
(V9)V/r) + ri + r2. 


(Pi^Ti + P 2 ^P 2 ) + div J' 


+ ^"^2^ m 


Thus, we identify JT” = -^^P^m{p)Vp-, and the equation for p becomes 

dtp + div {pv) = div (j^^^m{p)Vp'j + Ti + r 2 . (2.28) 

Remark 2.1 (Reformulations of the pressure and Darcy’s law). In the above derivation, 
we may consider the following pressure-type functions: 

• Let q-=p- A'I>{p) - P so that X-u = q + e- pp - N^^-a and 

V = K{V{-q-§ \Vpf) - BApVp) = -K{Vq + B div {Vp ® Vp)). (2.29) 


• Let p '■= p + N{p, a) so that X^ = p - pp - and 

V = K {V{N{p,a) + A'^{p) -p) - BApVp) 

= -K{Vp- pVp- N^^Va). (2.30) 
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(2.31) 


• Let p-=p + N(ip, a) - pp - so that Xv = p and 

V = K{V{N{(p, a) + A'^{p) - pp - N^aCr -p) - BApVp) 
= -K{Vp + pV p + aVN^a). 


We point out that (2.29) can also he obtained from the momentum balance of the Navier- 


Stokes-Cahn-Hilliard equations 

dt{pv) + div (pv 0 1 ))- div (p (Vv + (Vt))^)) + Vg ■ 


•div {BVp 0 Vp) 


by neglecting the inertia terms and replacing the viscous term with a multiple of the velocity. 
This is consistent with the classical derivation of Darcy’s law. 


Meanwhile, in (2.30) we have the gradient of the primary variables {(p,cr) multiplied by 


their corresponding chemical potentials {p, N^„), and vice versa in (2.31) (compare with the 


interfacial term K in Section 3 of and Eg. (2.34) of f^). It is common to reformulate 
the pressure as above to obtain equations of momentum balance in the Navier-Stokes- 
Cahn-Hilliard equations or the Cahn~Hilliard~Darcy equations that are more amenable to 
further analysis. See for instance 


2.4 Specific models 

2.4.1 Zero excess of total mass 


Assuming r 2 = -Fi =: T, so that there is no source term in (2.28), and let 


so that 


1 1 _ 1 1 
PS=— + —, 
P2 Pi P2 Pi 


Then (2.25) becomes 


r„ = ar, r,^ = p5r. 


divr) = oT, 

V = -K{Vp - pS/ip + N^^Vip), 
dtp + div (vp) = div (m(p)Vp) + psE, 
p = A4''(y)) - BAp + 
dt(T + div {(Tv) = div {n{p)\/N^fj) - S. 


(2.32) 


(2.33a) 

(2.33b) 

(2.33c) 

(2.33d) 

(2.33e) 


In the case that the densities are equal, i.e., pi = P 2 = p, then, a = 0 and ps = ^, 
becomes 


and (2.33) 


divi) = 0, 

(2.34a) 

V = -K{Vp - pVp + N^^Vp), 

(2.34b) 

dtp + V ■ Vp = div {m{p)Vp) + |r. 

(2.34c) 

p = A^' (p) - BAp + 

(2.34d) 

dta + V ■ VcT = div {n{p)VN^cr) - S. 

(2.34e) 
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2.4.2 Absence of nutrients 


Setting a = N{a,(p) = 0, then (2.25) simplifies to 


divt) = Pi^Ti + P 2 ^^ 2 -, 

V = -K{Vp - pV(p), 

dt(p+ div(r7(/9) = div{m{(p)Vp) + P 2 ^T 2 
p = A^'{(p) - BA(p. 


(2.35a) 

(2.35b) 

(2.35c) 

(2.35d) 


2.4.3 Zero velocity, zero excess of total mass and equal densities 

Suppose the volume-averaged mixture velocity v is zero, the excess of total mass Fi + r 2 
is zero and the densities are equal. Then, substituting d = 0 in (2.34) and neglecting the 
Darcy system (2.34^,b), we obtain 


dtip = div {m{(p)Vp) + , (2.36a) 

p = A^'((f) - BAip + (2.36b) 

dtcr = div {n{ip)VN^(j) - S. (2.36c) 


2.4.4 Boundary conditions for velocity and nutrient 

For the nutrient, we may prescribe a Robin type boundary condition: 

{n{(p)VN^a) • = c((Too - O') on dD, (2.37) 

where c > 0 is a constant, and o^o denotes a given supply at the boundary. When c = 0, 
we obtain the zero flux boundary condition: 

■ v> = 0 on dQ. (2.38) 

If we formally send c ^ oo, then we obtain the Dirichlet boundary condition: 

<T = cJoo on dD. (2.39) 

We may consider a boundary condition for the normal component of the velocity (which 
corresponds to a Neumann boundary condition for the pressure): 


-V ■ u = KVp ■ V = g2 on dD, (2.40) 

for some given function g 2 . We point out that a compatibility condition is required to 
hold if we consider the boundary condition (2.40) for the Models (2.25), (2.33), (2.34), 
and (2.35). Namely, if the mass exchange terms Fi and F 2 are given, then we require that 
<72 satisfies 


[ 92 ^ = f V u dH'^ ^ = [ div v dx 

Jan Jan Jn 


J pj^^Fi + 792^^2 dx for Models (2.25), (2.35), 


f aF dx 
Jn 


for Model (2.33), 


for Model (2.34). 
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However, if the mass source terms Fj depend on a or n, then considering (2.40) as a 
boundary condition would imply that (p, p and a have to satisfy 


f +P2^r2(<y9,(T,/r)dx= f -^2 d?^'^ ^ 

JQ J a fl 

Alternatively, we can prescribe a boundary condition for the pressure. Recall the refor¬ 
mulated pressure p and the Darcy’s law (2.30). We can prescribe a Dirichlet boundary 
condition; 


p = <71 on dQ, (2-41) 

for some given function gi. We may also consider the mixed boundary condition as in 
Section 2.3.3 of [12] (which corresponds to a Robin boundary condition for the pressure): 

ap-hv ■ u = ap + bKS/p ■ v - hKN^^Va • i/ = (73 on clD, (2.42) 


for constants a, 6 > 0 and a given function g^. 


2.5 Comparison to other models in the literature 
2.5.1 Absence of nutrients 

Scaling mass, permeability, and mobility appropriately, by setting 


ri = o, r:=r2, p2 = pi = i., k=i, m{ip) = i 


in (2.35), we obtain the following system 


divi) = r, 

V = -Vp + pVp, 
dtp + div (vp) = Ap + r, 

p = A^'i^p) - BAp. 


(2.43a) 

(2.43b) 

(2.43c) 

(2.43d) 


The existence of strong solutions in 2D and 3D have been studied in m for the case 
r = 0. For the case where F ^ 0 is prescribed, existence of global weak solutions and 
unique local strong solutions in both 2D and 3D can be found in [32]. We also refer the 
reader to [9] for the study of weak solutions to a related system, denoted as the Cahn- 
Hilliard-Brinkman system, where an additional viscosity term gdivD(v) is added to the 
left hand side of the velocity equation ( 2 . 43 | 3 ) and the mass exchange F is set to zero. 
Here, D{v) = |(Vr’ + (S/vY) is the rate of deformation tensor and rj is the viscosity. 


2.5.2 Zero velocity, zero excess of total mass and equal densities 

We consider the model (2.36) with the rescaled density p = I- Let V, A, C, Xa, Xv> be non¬ 
negative constants. For physically relevant values of the model variables, i.e., p e [-1,1] 
and cr > 0 , we choose 


r = {Va-A)h{p), (2.44a) 

N{p,(r) = + Xv^(^Y-p), (2.44b) 

S = Cah{p), (2.44c) 

where h{p) is an interpolation function with h{-l) = 0 and h{l) = 1 . 
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We have elaborated on the physical motivations for the particular forms of T and S in 
Section For the choice of N{ip,a), if both Xyp and Xa are positive constants, then for 
physically relevant parameter values, i.e., u > 0, and (p e [-1,1], 

= Xacr + XvC^ - ^ 0. (2.45) 

Thus, this choice of the flux VA^^o- provides two transport mechanisms for the nutrient a. 
The first term Xo-Vcr results in a diffusion process along negative gradients of a, while the 
second term is a chemotactic term that drives the nutrient towards the tumour cell 

regions. In particular, in the tumour cell regions {ip = +1}, the nutrient will only experience 
diffusion, while in the healthy cell regions {p = -1}, the nutrient will experience diffusion 
and active transport to the tumour. 

We point out that for this particular form of together with the zero Neumann 
boundary condition for p, we have 


ViV^o- ■ V = XcrVcr • V - ^ = Xo-Vcr ■ t> on dVl. 

With these choices, ( 2.36| ) becomes 

dtp = diY {m{p)Vp) + 2{Va - A)h{p), 
p = A^’{p) - BAp - x^cr, 
dta = div {n{p){xa'^cr-Xifi'^v)) -Ccrh{p). 


(2.46a) 

(2.46b) 

(2.46c) 


We remark that (2.46) is similar to Eq. (68)-(73) of [TT| . the two-phase diffuse interface 
tumour model (Eq. 5.27) of |38] . and model A42 of [28]. The only difference between these 
three models and ( 2.46| ) is that the flux for the nutrient equation (2.46c) consists of an 
advection term and a Fickian diffusion term for [T^, while in |38( 128] . the nutrient is in 
a quasi-steady state and the flux for the nutrient equation is a Fickian diffusive flux. We 
point out that in [TUEEIEHI, h{p) is replaced by p in the definition of F and S. Since, in 
their notation, p e [0,1] denotes the tumour volume fraction instead of the difference of 
volume fractions. 


Next, choosing N{p,a) as in (2.44b) above, and 


T=^P{p){N^-p), S = P{p){N^-p), 


where P{p) is a non-negative function, then ( |2.36 ) becomes 

dtp = div {m{p)S/p) + P{p){xaCr + x^O- - 
p = A^'{p) - BAp - x^cT, 

dta = div (n(<y9)(xaVcT - x>f^p)) - P{^){Xaa + Xv{^-v) - p)- 


(2.47a) 

(2.47b) 

(2.47c) 


This is similar to the model derived in |29j , where the chemical potentials and p enter 


as source terms in (2.47a) and (2.47c). The specific form for F is motivated by linear 
phenomenological constitutive laws for chemical reactions. The non-negative function 
P{p) takes on the form 


P{^) 


[(5Po( 1 + 7’), if </?>-!, 


lO, otherwise, 

for positive constants 5 and Pq. Subsequently, if we choose 

X(T = 1, Xv = = 1 


(2.48) 
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in (2.47), we obtain 


dt‘f = A^ + P{ip){a - ij), 
jj. = A^'{(p) - BAip, 
dta = Act - P{(p){o- - n). 


(2.49a) 

(2.49b) 

(2.49c) 


This is the model studied in |22] . for a more general function P{^p) than (2.48), while a 
viscosity regularised version of (2.49) (where there is an extra adtfi term on the left hand 
side of ( 2.49a| ) and an extra adtip term on the right hand side of (2.49b) for a positive 
constant a) is studied in m- A formal asymptotic limit for the viscosity regularised 


version of (2.49) is derived in m- 


3 Sharp Interface Asymptotics 


We consider Model (2.25) with the following choices and assumptions: 

Assumption 3.1. 

• A = ^ and B = fde for positive constants j3,e > 0. 


• N((p,a) is chosen as in (2.446) with constant positive parameters > 0- 

• The mass exchange terms Tj, i = 1,2, and the nutrient consumption term S depend 
only on a, p, and ip, and not on any derivatives. 


• The mobilities m{(p) and n(ip) are strictly positive and continuously differentiable. 

• The potential ih is chosen to be either the smooth double-well potential '!'((/?) = j(l- 
<y9^)^ or the double-obstacle potential 




< 1 , 


if 

otherwise 


(3.1) 


With these choices. Model (2.25) becomes 


divu = p^^Ti{a,ip,p) + p2^T2{a,ip,p), 

V = -K{Vp - pVip - x<pCrV(/7), 

dtip+ div{vip) = div(m((/?)V/u) +p2^T2{a,ip,p) -pl^Ti{a,ip,p), 
p = ^d>'{ip) - feAip - Xpcr, 

dtcr+ div((Tr>) = div (n(<y 9 )(XffVcj - x^^Vv^)) -S{a,ip,p). 


(3.2a) 

(3.2b) 

(3.2c) 

(3.2d) 

(3.2e) 


We point out that in the case of the double-obstacle potential, the “derivative” T' is to 
be understood in the sense of subdifferentials, i.e.. 


'^'{ip) = -ip + 9/[_i^i](v?), 5/[_i_i]((/?) = 


'(-oo,0] iiip = -l, 

0 if \ip\ < 1, 

[0, -i-oo) \i ip = +1, 


(3.3) 


and ( 3.2d[ ) will have to be formulated in terms of the following variational inequality: 

r 0 

/ -/r(V'-¥?)--</?(V’-V^) + /5eV¥5-v('0-(/::)-X¥^cr(V’-</?)dx>O, (3.4) 

JQ £ 
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for all Tp ^K, '■= {r] e : \r]\ < 1}. 


We perform a formal asymptotic analysis on Model (3.2) in the limit e ->• 0. Details of 


the method can also be found in [3 [3 miss]- The following assumptions are considered: 

Assumption 3.2. 

• We assume that for small e, the domain D can be divided into two open subdomains 
D=‘'(e), separated by an interface S(e) that does not intersect with dQ. 


• We assume that there is a family {(ps,v^,p^, yLeT(ye)e>o of solutions to (3.2), which are 
sufficiently smooth and have an asymptotic expansion in e in the bulk regions away 
from S(e) (the outer expansion), and another expansion in the interfacial region 
close to S(e) (the inner expansion). 

• We assume that the zero level sets of ips converge to a limiting hypersurface Sq 
moving with normal velocity V. 

The idea of the method is to plug the outer and inner expansions in the model equations 
and solve them order by order, in addition we have to define a suitable region where these 
expansions should match up. 


We will use the following notation; (3.2e)Q and (3.2e)“ denote the terms resulting 
from the order a outer and inner expansions of (3.2e), respectively. 


3.1 Outer expansion 

We assume that for e the following onter expansions hold: 


fe - fo + ^fi + +- 


To leading order (3.2^)^ gives 


-/3^'ipo) = 0 . 


(3.5) 


The solutions to (3.5) corresponding to minima of T are po = ±1, and thus, we can define 
the tumour tissues and the healthy tissues region by 

Qt := {x e D : (po{x) = 1}, := {x e D : (po{x) = -1}. (3.6) 


Then, thanks to Vpo = 0, we obtain from the equations to zeroth order: 


divuo = p^^Ti{ao,ipo,po) + P2^'['2{(ro,(po,po), (3.7) 

Vo = -KVpo, (3.8) 

- div {m{po)Vpo) = P 2 ^(l - Vo)^ 2 {cro,Po, Po) - ^1^(1 + y?o)ri((To, (/?o, /^o), (3.9) 

9tCJo + div((ToUo) = div(n((^o)Xf7Vcro)-5 ((To,(/::o,/^o)- (3.10) 


For the double-obstacle potential, we obtain from (3.4)^^, 


y' -/3(^o(V’o “ T^o) dx > 0 for all V’o ^ ^ ^ F7^(D) : \p\ < 1}. 


For this to hold for all |'0o| ^ 1; we require that po - ±1; a-nd thus we can define and 
VLh as before, and also recover (3.7), (3.8), (3.9), and (3.10). 
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3.2 Inner expansions and matching conditions 

By assumption, Sq is the limiting hypersurface of the zero level sets of In order to 
study the limiting behaviour in these parts of we introduce a new coordinate system. 

We introduce the signed distance function d{x) to Eg, and set z = ^ as the rescaled 
distance variable, and use the convention that d{x) < 0 in Qh, and d{x) > 0 in Thus, 
the gradient Vd points from Qh to Qt, and we may use Vd on Eq to denote the unit 
normal of Eq, pointing from to 

Let g{t, s) denote a parametrization of Eq by arc-length s, and let iz denote the unit 
normal of Eq, pointing into the tumour region. Then, in a tubular neighbourhood of Eq, 
for sufficiently smooth function /(x), we have 

f{x) = s) + ezu{g{t, s))) =: F{t, s, z). 

In this new (t, s, z)-coordinate system, the following change of variables apply, compare 

IS!: 

dtf =--VdzF + h.o.t., 
e 

Vx/= +VsoT+ h.o.t., 

where V is the normal velocity of Eq, Veo 5 denotes the surface gradient of g on Eq and 
h.o.t. denotes higher order terms with respect to e. In particular, we have 

A/ = diVa;(Vx/) = \dzzF+ - divj^g{dzFu) + h.o.t., 

= -Kd^F 

where k = - div is the mean curvature of Eq. Moreover, if t; is a vector-valued function 
with V(t, s, z) = v(x) for x in a tubular neighbourhood of Eq, then we obtain 

diva,t; = ^dzV ■ v -i- divEoI^ + h.o.t.. 

We denote the variables ipe, Ve, Pe, Ps, (Xe in the new coordinate system hy Ve, Pe, 

Ce, respectively. We further assume that they have the following inner expansions: 

F^{s,z) = Fo{s,z) + £Fi{s,z) + ..., 

for F^ € {<I>£, Pe, He, Ce}. The assumption that the zero level sets of converge to Eq 
implies that 


$o(Ls,^ = 0) = 0. (3.11) 

Furthermore, we make the following assumption: 

Assumption 3.3. For the double-obstacle potential, we assume that the inner variable $£ 
is monotone increasing with z and the interfacial layer has finite thickness of 21, where the 
value of I will be specified later. For the double-well potential, we take 1 = oo. Furthermore, 
we assume that 


S,Z = +1) = +1, S, Z = -1) = -1. 


(3.12) 
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In order to match the inner expansions valid in the interfacial region to the outer 
expansions of Section 3.1 we employ the matching conditions, see 


lim Fo(t,s,z) = fo{t,x), 
z^±l 

lim dzFQ{t, s, z) = 0 , 

lim dzFi{t,s,z) = V/q (t,x) ■ ly, 
z^±l 


(3.13) 

(3.14) 

(3.15) 


where /q (t, x) ■= lim^^o /o(i) x ± 5u) for x e Sq. Moreover, we use the following notation: 
Let (5 > 0 and for x e Sq with x - ^ and x + 6u ^ we denote the jump of a 

quantity / across the interface by 

[f]Jj '■= lim/(t, X + 5v) - lim/(t, x - 6u). (3.16) 

( 5\0 ( 5\0 


For convenience, we define the constant 7 > 0 to be 

1 , A , . . 2\/2 


7 


/ oo ^ 

-sech^(z/V 2 ) dz = —— for the double-well potential, 

00 2 3 

J " cos^iz) dz = ^ for the double-obstacle potential. 


(3.17) 


3.2.1 Expansions to leading order 

To leading order {3.2d)]^ gives 


dzz^O - ^'($ 0 ) = 0 . 


(3.18) 


Using (|3.11|) we obtain that $0 can be chosen to be independent of s and t, i.e., <I>o is only 

(3.19) 


a function of z, and solves 

4>(;(z)-T'($o(^)) = 0, ^o(0) = 0, 4>o(±/) = ±l. 

For the double-well potential, we have the unique solution 


<l>o(z) = tanh I 

\^/2 


(3.20) 


Furthermore, multiplying (3.19) by ‘hg(z), integrating and applying the matching condi¬ 
tions (3.13) and (3.14) to 4>o gives the so-called equipartition of energy: 


- |‘hQ(z)p = 4'(<I>o(^:)) for all \z\ < 00 . 




Similarly, for the double-obstacle potential, we obtain from (3.4)} , 

< 1 . 


-/3($o + dzz^o){ilJ - $ 0 ) dx > 0 for all 


(3.21) 


(3.22) 


For (3.22) to be satisfied, it suffices to consider <l>o as a function only in z which solves 
$o(^) + $o(2) = 0, $o(0) = 0, $o(±0 = ±l- (3-23) 
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A solution to (3.23) is 




+ 1 

sin( 2 ;) if \z\ < 

-1 if 2 ;<-|, 


(3.24) 


so that / = I for the double-obstacle potential, and we deduce from (3.12) that for the 
double-obstacle potential, 


s, ±5) - 0 . 


(3.25) 


Moreover, we obtain the equipartition of energy (3.21) via a similar argument to the 
double-well potential. Thanks to the equipartition of energy (3.21), and the definition of 
7 (3.17), we point out that 

dz = 2 T( 4 >o( 2 ;)) dz = 7 . (3.26) 

For the rest of this section, we do not differentiate between the two cases of potentials, 
and use the notation that I = | represents the case of the double-obstacle potential and 
I = 00 repr esents the case of the double-well potential. 

Next, (3.2a)y^ gives 


dzVo ■ 1 / = 0 . 


(3.27) 


Integrating from -I to I with respect to z, and applying the matching condition (3.13) to 
Vq yields 


We have from {3.2c)j^, 


[vo]h ■iy-=VQ-u-VQ-u = 0 . 


dzim{<^o)dz'^o) = 0 . 


(3.28) 


(3.29) 


Upon integrating and using the matching condition (3.14) applied to Hq, we obtain 

m{^o)dz^o{i, s,z) = 0 for all \z\ < 1. 

Since |4>o(2:)| < 1 for \z\ < I and m(<I>o) > 0, we have 

dz'B(j{t,s,z) = 0 for all \z\ < 1. (3.30) 

Thus, integrating once more with respect to z from -I to I, and applying the matching 
condition (3.13) to Hq, we obtain 


[i^o]h - d - 


To leading order, the nutrient equation (3.2e)j yields 

dz{n{^o)x<TdzCo) - {n{^o)Xip^o{z)y = 0 . 


(3.31) 


(3.32) 


Integrating and using the matching condition (3.14) applied to both Cq and <I>o leads to 


n{^o){XadzCo - 4 * 0 ( 2 :)) = 0 for all | 2 ;| < 1. 
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As n($o) > 0, we see that 

XcjdzCQ{t,s,z) = Xip^'oiz) for all jz] < 1. (3.33) 

Integrating once more with respect to z from -I to I, and applying the matching condition 
(3.13) to Co and $o then gives 


Lastly, {3.2b )yields 


X(7 Xa 


QzPq = (So + X(/jCo)<ho. 


Integrating and applying the matching condition (3.13) to Pq and Ho leads to 

T 

[po]h = + Xv C'o(L s, z)^q{z) dz. 

Thanks to ( 3.33[ ), we see that 

/''co<I>;dz = ^ CcodzCodz = ^ f'dzl^-^]dz 
J I x<p ^ X^p ^ \ ^ / 


(3.34) 


(3.35) 


(3.36) 


X. 


Then, (3.36) becomes 


, [ICoP]', —NAh' 


[poIh = 2/^0 + ^ [kol^]^ • 


(3.37) 


(3.38) 


3.2.2 Expansions to first order 


For the double-well potential, to first order, we obtain from (3.2(7)^, 


(3.39) 


, 0 T"($o)^i - fidzz^i + - XpCo = Ho. 

We multiply ( 3.39[ ) with <I>g and integrate with respect to ^ from -oo to oo, which gives 

/ oo 

'P‘o{t, s)‘^q{z) dz 

OO 

/ oo n 

/3(T'($o))'^i-/35,.$i4>(, + /3k|$[,| -XpCo^odz. (3.40) 

oo 


Applying integration by parts and the matching conditions (3.13) and (3.14) applied to 
‘ho, and using that 'I''(±l) = 0, we see that 

X oo 

(T'($o))'^i -5^^^i$odz 

oo 

X oo 

a^$i(4''($o)-$o) 

__ oo V 


dz, 


=0 by ( |3.13| |,( |3.14| | 


=0 by ( |3.19| | 


and so the first two terms on the right hand side of (3.40) are zero. Then, using (3.26), 
(3.30), and (3.37), we obtain from ( |3.40[ ), 


2^0 = I3kx-Xp f C'o$odz = /37K-^[|cro|^]^. 

..y “■ oo ^ 


(3.41) 
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Moreover, together with (3.38), we obtain 


[po]h = /37«^- 


(3.42) 


Meanwhile, for the double-obstacle potential, to first order, we obtain from (3.4) 


/! 


^(-So - XypC'o - Pdzz^i - - 4>o) dx > 0 for all \^p\ < 1. (3.43) 


Since |<l>o(.2;)| < 1 for \z\ < we can test with -0 = 4>o + A with either non-positive or 
non-negative A e /C, leading to the equality 


-Ho - Xv>Co - f3dzz^i - (3^1 + nf3^o = 0. 


Multiplying with 4 >q and integrating with respect to 2 ; from - ^ to ^, and applying matching 
conditions leads to 


r' — r' — r' — 

-2fio + l3K j I dx-j ^ x<fiCo^odz = P j ^ dzz^i^o +^i^'odz. (3.44) 


Upon integrating by parts and using (3.23), the matching conditions (3.14) for <f>o, (3.15) 
for $ 1 , and (3.25), we see that 

dzz^i^o + $i^>o dz = - j_l (4>(; + dz = 0. 


Then, using (3.26), and (3.37), we obtain from (3.44) the following solvability condition 
for $ 1 : 


2/ro = /37«^- Y [koP]^- 

Lastly, thanks to (3.30), we obtain from (3.2c)k and ( |3.2e )k) respectively, 

(-V + Vo ■ iz) 4 >o = dz{m{^o)dz^i), 


(3.45) 


and 


(-V + Fo • r^)dzCo 

= dz{n{<^o){XadzCi - Xipdz^i)) 

+ (9^(n'(^>o)^i {XadzCo - X^^'o)) + divso(?^(^o) {XadzCo - X^^K) 

=0 by ( |3.33D =0 by ( |3.33| | 

= dz{n{^o){xadzCi - Xipdz^i))- 


(3.46) 


Thanks to (3.27), upon integrating from -I to I with respect to z, and applying the 
matching condition (3.15) to Hi, we obtain from (3.45) 

2(-V + no-iz) = [m(v?o)V/Uo]^-iz. (3.47) 

Similarly, thanks to V(/9o = 0, upon integrating from -I to I with respect to z, and applying 
the matching condition (3.15) to Ci and 4>i, we obtain from (3.46) 

(-V + vo-v) [ctq]^ = x<7 [n{(po)Vcro]H ' (3-48) 
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3.2.3 Sharp interface limit for Model (3.2) 


In summary, we obtain the following sharp interface limit from Model (3.2): 


Vo = -KVpo 

in rij' u XIh , 

(3.49a) 

divr7o = /5i ri((To , 1,/To ) ■*■ 6*2 r 2 (o'o ! /^o ) 

in ri'T, 

(3.49b) 

divno = Pl^ri(fj^,-1, ^^) -H P 2 ^r 2 ( 0 -^,-1,/U^) 

in Qh, 

(3.49c) 

-m{l)Apl = 

in XIj’, 

(3.49d) 

-m{-l)Ap^ = 2p2^T2{a^ ,-l,Po) 

in Qh, 

(3.49e) 

dtal -h div(cj^no) = n(l)x^ACT^ - 5((t^, 1,/r^) 

in XIt, 

(3.49f) 

dtao + div {a^vo) = n(-l)xaA(T^ -S{ao, -1, Po) 

in Qh, 

(3.49g) 

together with the free boundary conditions 



[no]^-i^ = 0, [cJo]^ = 2^, [po]h = ^7'« 

X(7 

on So, 

(3.50a) 

= 2/io =/37 k-y [kop]^ 

on So, 

(3.50b) 

2(-V + vo-u) = (m(l)V/Uo - m{-l)Vpo ) ■ 

on So, 

(3.50c) 

2—(-V + V 0 -V') = Xa{n{l)\/ao - n(-l)Va^) • u 

Xa 

on So, 

(3.50d) 

where 7 is defined in (|3.17D. Note that we can write 



[koP]^ = Mh i^o + (^ 0 ) =■ 2^0 [c^oIh > 


(3.51) 


where (Jq := denotes the average of the nutrient concentrations from both sides 


of Sq. Thus, using (3.34), we can rewrite (3.41) and (3.506)2 as 


1^0 = ]^hi^-croX^ on So. 


(3.52) 


3.3 Specific sharp interface models 

In this section, we take T as the double-well potential. 


3.3.1 Sharp interface limit of the new active transport model 

Choosing as before N(ip,a) = ^ \af' + X(^cr(l - ip) and 


r(cr, p) = {Va - A)h{p), S{a, p) = Cah{p), 
m(p) = mo > 0, n{p) = no > 0, 


(3.53) 

(3.54) 


for some positive constants V, A and C in Model ( |2.33 ), we obtain the Cahn-Hilliard- 
Darcy model (1.1) in Section]^ Then, the sharp interface limit of Model (2.33) with (3.53) 
is given by 


-Apo = 
-moApo = 


|o in Q.H, 

Ups-a){Vao-A) in fir, 
|0 in hl/f, 


(3.55a) 

(3.55b) 
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dtCTo - div (KaoVpo) - noX^Acro = 

nT n 


1 -CcJo in XIt, 
jo in Qh, 

(3.55c) 

[cro]^ = 2—, [po]^ =/37 k on So, 

(3.55d) 

Xa 


Mo + ^oX<fi = on So, 

(3.55e) 

T 

"lo [^Mo]h ■ ^ ^0, 

(3.55f) 

T 

noXa [V(To]j|^ ■ V on Sq. 

(3.55g) 


-2^{V + KVpo-i^ 

X-7 

The active transport term n{(p)'\/{xip^) in the flux for the nutrient results in the jump 
term 2— in (3.55d)2 which is a new feature of the proposed model. 

Xcr 1 I 

3.3.2 Linear constitutive laws for chemical reactions 


Let us consider Model (2.47) with m{ip) = n{ip) = 1, and Pilf) defined as in (2.48). Then, 
we obtain that 


P{^o){Xa(ro + X</j(l - <^o) - Ato) = 


1 2 ( 5 Po(Xf 7 ' 7 o - Mo) in fir, 
jo in Qff. 


(3.56) 


This was the setting introduced in [29]. Hence, we obtain from (2.47) the following sharp 
interface limit: 


-Apo = 1 

|25Tb(XaO-o -Po), 

in Qt, 
in Qh, 

(3.57a) 

dtCJO - XfrAcTo = 1 

\-26Po{xaCro- Mo), 
[0 

in Qt, 
in Qh, 

(3.57b) 

\ T T 

Mo = - ^0X<fi, [mo]h = 0) [(^o]h = 2— on Sq, 

^ Xct 


(3.57c) 

-2V = [V/ioln ■ -2—V = Xa [VcTo] J on Sq. 


(3.57d) 


We point out that the diffuse interface model studied in m takes a different mass transi¬ 
tion term T and a different consumption term S. More precisely, the following choices are 
considered: 


r = ^Pi^)icr -5p), S= ^P{ip){a - 6p), 


where 0 


P(^) := 


\2Po^{(p) if [-1,1], 
0 otherwise . 


*In [30], the choice of T((p) is actually 


Pi^) := 


|2Po(T((p))i if [-1,1], 
[O otherwise . 


(3.58) 


(3.60) 


(3.59) 


This presents some difficulties in the analys is of t he out er expansions, as P'(±l) t 0 for the choice = 
{1 - <p^)^. Thus, we do not recover (3.61a| and (3.61bl. However, the formal analysis in [30| is different 
compared to what we present here, and it turns out that the analysis in m has to be modified and will 
only work for (3.601. 
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The term acts as a regularisation on the HausdorfF measure restricted to the limiting 

hypersurface Sq, and hence the reaction term a-6^ will appear in the interfacial relations 
of fjL and a rather than in the bulk equations. More precisely, we obtain 


A|Uo = 0 in fly u (3.61a) 

dtCTo = A(To in fl'r u (3.61b) 

= 0, [fJo]^ = 0, 2^0 = iS'yn on So, (3.61c) 

-2V = [V/io]^ • + Po'y{(^o - SiJ-o) on So, (3.61d) 

0 = [Vo-o]^ • ly - Po^ido - 5/ro) on So, (3.61e) 



3.3.3 The limit of vanishing active transport 


We con sider Model (2.46) with a quasi-steady nutrient (i.e., neglecting the left hand side of 
(2.46c)), with positive constants D and A, and the interpolation function h{ip) = |(1 + (p), 
we set 


V{^) ^ = 1(1 + I)) + |(i _ z)), (3.62a) 

mi^p) = ^{1 + (p)^, n{ip) = XV{ip)x~^= (3.62b) 

so that, we obtain 

dtp> = d\Y {^{1 + ip)‘^Vp)+Va{(p + 1) - A{p} + 1), (3.63a) 

/i = ^T'(z?)-/3eA(^-X<^o-, (3.63b) 

0 = div {'D{ip)S/cr) - A div {'D{ip)S/(p) - -Ca{ip + 1). (3.63c) 


The specific choice (3.62) allows us to control the influence of the active transport term 
n{ip)xtpX7(p via the parameter A, while preserving the chemotaxis term -X(fid in (2.46b). 
Hence, we have “decoupled” chemotaxis and active transport. 

Moreover, if we consider Eq. (68)-(70) of [H] with the choice cj)= |(1 + y?), ^ = 1, and 
the rescaling p. ep, the resulting phase field model almost coincides with Model (3.63) 
with the exception of the additional term X div ('D(ip)\/ip) in (3.63c). 


We briefly state the derivation of the sharp interface limit for Model (3.63). From 
(3.635 )q^ we have (po = ±1 and the domains Qt and ^h- From (3.630)^ and (3.63c )q we 
obtain 


0 = div(|(l + (po)‘^Vpo) +Vaoiipo + 1) - A{ipo + 1) in XIt u 
0 = div (2 ?((/?o)V<to) - ^Cao((pQ + 1) in 12^ u 


From the leading order inner expansion (3.635 )^^, w e obtain (3.18), and subsequently the 
profile (3.20) and the equipartition of energy (3.21). From (3.63c)J^ we have 


d,{V{^o)d,Co - AP($o)$o) = 0. 
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Integrating and nsing the matching conditions (3.14), we obtain 


V{<l>o){d,Co-X<^o) = 0- 
Since P(<I>o) > 0 for |4>o| < 1, we obtain 

dzCo{t,s,z) = A<I>q( 2 ) for all \z\ < 
and upon matching, we obtain 

dz{{l + ^ofdz'^Q) = Q. 


(3.64) 


(3.65) 


While from (3.63a we obtain 


Integrating and using the matching condition (3.14) applied to Hq we deduce that 

(1 + <ho('2;))^^z^o(i) s, z) = 0 for all \z\ < oo. 


Since |4>o(2:)| < 1 for \z\ < oo, we obtain that dz'^oit, s, z) = 0 for \z\ < oo. 

)?, 


To first order, we obtain from (3.636)5, 


Ho = /34''($o) - Pdzz<^i + f3K<^o - X^Cq. 

Multiplying by ‘hg, using that Hq is independent of and applying integration by parts 
and matching conditions (3.13[) and (|3.14) to $ 0 ; we obtain, in the same spirit as (3.40), 


2^0 = dz = ^ ^ [koP] 


,2W 
iH 


where we have used (3.64). Applying (3.51), and (3.65), we see that 

2/io = /37 k - [cTo]^ = (3-fK - 2x^ao, 

where we recall that ao ■= + (7^) is the average of the nutrient concentration at the 

interface. Meanwhile, thanks to (3.64) we obtain from (3.63c )k, 

0 = dz{V{^Q){dzCi - A9,$i) + P'(^o)^i(9,Co - X^'o)) 

= dziv{<^o){dzCi-xdz<^i)). 


Integrating with respect to z from -oo to oo and applying the matching condition (3.15) 
to Cl and <hi leads to 

0= [T>((/?o)Vcro]^ ■ iz. 


Lastly, thanks to = 0, we obtain from (3.63a)j^, 

-V^’o^^dz{{l + <^ofdzEi). 


Integrating from -oo to oo with respect to z, and applying the matching condition (3.15) 
to Hi gives 

-2V = 2V/ao • u. 
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Thus, the sharp interface limit of Model (|3.63) is 

Val- 
\Cao 




A(T( 


0 


0 


A in VLti 
in Qy, 
in Q.H, 


[o-o]^ = 2A, 2/xo =/37 k-Xvp(o-^ + o-?) on So, 
0 = (V(J(f-T>Va^)-i^, -V = V/Xo ■ ^ on So. 

In addition, we can use (3.66c)i to rewrite ( |3.66c )2 as 

2/io = - X<pi‘^cro - 2A). 

Next, sending A ^ 0 in (3.66) leads to 

[o-o]// = 0, 2^io = [3 jk-2x^(To, 
and we define the bulk velocity and pressure via the relations; 

-V(p-X<^cro), P ■= ^^o + Xv^cro. 

Then, we deduce that 

V = -V^o ) div V = -A/Tq = - Al) in 0^, 


and from ( |3.66d[ ) and ( |3.68[ ), 

-V = V/Uq v = -vu = V(p- Xv^o-q) ■ ^ on So, 
P = P0 + XipCTQ 


1 


Thus, we obtain 


-/37k on So. 


div V = Vuq - M in VLt, 

V = -V(p- Xv^o-q) in AIt, 
I Cao in Xlx, 

lO in A2 h, 


Ado 




(VcrJ - DVa^) • iv = 0 on So, 


P- 


^/37k on So, 


-Vp • V + X(/jVdo • IV = V on Sc 


(3.66a) 

(3.66b) 

(3.66c) 

(3.66d) 

(3.67) 

(3.68) 

(3.69) 

(3.70) 

(3.71) 

(3.72) 

(3.73a) 

(3.73b) 

(3.73c) 

(3.73d) 

(3.73e) 

(3.73f) 


which coincides with the sharp interface model (Eq. (79)-(81), (83)-(86)) of [H]. We point 
out that the same sharp interface limit (3.73) can be recovered if we set A = e in (3.63). 
We introduce the parameter A in (3.63) in order to study the effect of active transport on 
the linear stability of radial solutions to (3.66), see Section 0 below. 

Let us also remark that the mobility m{p) = |(1 + P) is degenerate in the region 
{if = -1}, and thus the bulk equation for in XIh and the interfacial condition for V-v 
on So remain undetermined in (3.66). Furthermore, if m(if) is chosen to be degenerate in 
the bulk regions {f = ±1}, then we obtain from the outer expansion (3.630)^ the following 
equations 

0 = Vaoifo + 1) ~ Al((po + 1) in u 

In particular, we see that (T((" = 4 is a constant in which is inconsistent with (3.63c 


o- 


Hence, it is necessary that the mobility m{f) is not degenerate in the tumour region 
{ip = l}. 
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4 Linear Stability Analysis 


Let us consider the sharp interface model (3.66). By sending the active transport param¬ 


eter A to zero, we recover the sharp interface model (Eq. (79)-(81), (83)-(86)) of [Hj. In 
this section, we extend the linear stability analysis of [HIS! to include the effects of active 
transport. For the linear stability analysis of a one-phase model, we refer to [MIES]. 

4.1 Radial solutions 


We now drop the index 0 in (3.66), and let = B/j(0) denote the d-dimensional ball, 
h = 2,3, of radius R centered at the origin. We assume that the interface S is a (d- 1)- 
sphere of radius q{t), partitioning the domain H into VLt and as follows: 


The outer unit normal ^{p) at a point p e S is given as 

( \ P P 

^KP) = r-i = 


while the normal velocity V is given as 


\P\ 


dt 


The mean curvature k for a (d- l)-sphere radius vq is given by 

d-1 


ro 


(4.1) 


(4.2) 


(4.3) 


where d denotes the dimension. Then, for radially symmetric solutions v?(|a;|) = 


/i(|a;|) = p(r), <t(|s|) = cr(r), (3.66) becomes 

d-1 , 


n , 


a" + 


r 

d-1 


= A- VaT in r < q{t), 

Ca in r < q{t)^ 
in r > q{t)^ 


- a = 

r 10 


[cr]^ = 2A, 2pT = h'Kw - XipicTT + cfh) on r = q{t), 


Qit) 


dq 


a'rp = Da'fj, = Pr on r = q{t). 

We complete ( |4.4[ ) with the following boundary conditions: 

anir = R,t) = CToo, arir = 0,t) < °°, PT{r = 0,t) < 


(4.4a) 

(4.4b) 

(4.4c) 

(4.4d) 

(4.5) 


where aoo denotes the concentration of a nutrient supply from the boundary dXl. 

Upon solving the differential equations and applying the interface and boundary con¬ 
ditions, we arrive at the following radial solutions: 


(JH{r,t) = 


[cTcx) + a 2 (t)(log(r) - log(i2)) ford =2, 


^ oo 


«3(i)(F-s) 


for d = 3, 


(4.6a) 
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aT{r,t) 


HT{r,t) 


h 2 it)hiAr) ford = 2 , 

r 


(4.6b) 

( A ^ V 

-r --b 2 it)IoiAr) + C 2 it) 

for d= 2, 

(4.6c) 

A 2 D, sinh(Ar) 

Mt) '+C3it) 

^ 6 C r 

for d = 3, 


where, for a e M, Ia{x) denote the modified Bessel functions of the first kind : 

oo ^ / r^\ 2k+a 


laix) = Y, 


k=ok\T{k + a + l) 


(ir 


(4.7) 


Here, r(') denotes the Gamma function. Together with the modified Bessel functions of the 
second kind, Ka{x) '■= f pairs {la^Ka} are the two linearly independent 

solutions to the modified Bessel’s equation: 


2 dy 2 2 

X —— + X— - X y = a y. 

dx dx 


(4.8) 


Moreover, for the case a = 0, the modified Bessel functions Io{x), Kq{x) satisfy the 
following properties 


d r 

/o(0) = 1, lim.Ko{x) = + 0 O, —Io{x) = Ii{x), / xlo(x) dx = x/i(x). 

dx J 


Furthermore, the coefficients in (4.6) are given as 
A^ = C, 
a2{t) = 


q{t)AIi{Aq{t)){a^ + 2A) 


DIo{Aq{t)) - Aq{t) \og{q{t)/R)Ii{Aq{t)) 

Rq{t){l - y(t)Acoth(Ag(t))) 


asit) = (cToo + 2A) 
62 (t) = 
hit) = 


(R - y(t))((?(t)Acoth(Ay(t)) - 1) + DR 
DicToo + 2A) 


DIo{Aq{t)) - q{t)A\og{q{t)/R)Ii{Aq{t)) 
((TtK, + 2A) DRqit) 


sinh(A( 7 (t)) {R - q{t)){q{t)Acoi\i{Aq{t)) - 1) + DR 


C2(t) = ^ + 




2 , Pi 


“6 "(*) %(*) 


+ (^ - Xv] ^>2(i)^o{A$(i)) , 
X \ 7 / X sinh(Aq(t)) 


and the differential equation satisfied by qit) is 


dt 


~Q+jb2it)hiAq) 

A P/Acosh(Ag) sinh(Ay)'^ 

--^9+ 3(G^I 


Q 


q^ 


I 


Qit) 

for d= 2, 
for d= 3. 


(4.9) 

(4.10a) 

(4.10b) 

(4.10c) 

(4.10d) 

(4.10e) 

(4.10f) 

(4.10g) 


(4.11) 


We point out that, thanks to the boundary condition ^^(r = 0,t) < 00 , the solution ax 
does not contain any terms involving K^iAr) (in d= 2) and cosh(Ar)/r (in d = 3). 
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4.2 Perturbation of radial solutions 

We now consider a perturbation of a radially symmetric tumour, whose radius w is given 
by 


w{r,e,(j),t) = q{t) + 6{t)z{e,4i), z{e,<i)) = 


[cos(Z 0 ) for d = 2 , 

\yi,m{0,(j)) ford = 3, 


(4.12) 


where q{t) is the radius of the unperturbed interface, 5{t) is a dimensionless perturbation 
size, is a spherical harmonic with I and 6 denoting the polar wavenumber and angle, 
and m and (j) denoting the azimuthal wavenumber and angle, respectively. We will denote 
the radial solutions in (4.6) by a^, and and consider 


aHir,e,(j),t) = + U{r,t)5{t)Z{e,(j)), 

arir, 6, </>, t) = t) + V (r, t)6{t)Z{e, 4>), 

+ W{r,t)6{t)Z{e,4>), 

where we assume that I^^t) solve ( 3.66| ). Therefore, we get 

A{iJ,^ + WSZ) = A-V{a^ + V6Z) in r<w, 

A{a^+ V6Z) = C{a^+ V6Z) in r<w, 

A{a*^ + U5Z) = 0 in r > rc, 

~ + {V - U)6Z = 2A on r = w, 

2(/iy + W5Z) = Pn'y - 2Xif{(yT + 2x^A on r = w, 

{(T^)r + 6^/{VZ) ■ If = D{{a*u)r + 5V{UZ) ■ u) on r = w, 


-^^-zf^={f^T)r + SV{WZ)-U 


on r = w. 


(4.13a) 

(4.13b) 

(4.13c) 


(4.14a) 

(4.14b) 

(4.14c) 

(4.14d) 

(4.14e) 

(4.14f) 

(4.14g) 


Here, we used the more convenient form (3.67) of (3. 660 ) 2 . Next, we linearise (4.14) 
about the original interface r = qto derive the equations satisfied by U, V, W and 6. We 
introduce the Laplace-Beltrami operator on the (d- l)-sphere, for d = 2,3: 


Cd 


502 


1 


52 


52 5 

502 ^ ^ do sin( 0)2 5(/>2 


for d = 2 , 
for d = 3, 


(4.15) 


so that the Laplace operator can be decomposed into 

= frr -I- fr + -^^df ■ 


(4.16) 


Moreover, the function Z{6,(l)) defined in (4.12) satisfies 


CdZ{9,ct>) = Ci4Z{0,(^), Ci,d 


{-1(1 + 1) 


for d = 2 , 
for d = 3. 


(4.17) 


From the bulk equation ( |4.14a ) we obtain 

A^iT + 5A(WZ) = A- Vut - 6VVZ, 
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and so, using that Afi^ = A - Va^, we deduce that 


VV + Wrr + -—-Wr + ^Ci,dW = 0 in r < g. 

For the interface conditions, we employ Taylor’s expansion and neglect terms of order 
0{5^). For instance, from (4.14d), we see that 

2A = (TT{q) + {(7T)r{q){w -q)- cr^iQ) " {^H)r{q){'w -q) + {V - U)SZ + 0(6'^). 

Then, by (3.66c)i, (3.66d)i, and ( 4.12| ), we obtain 

U{q,t) -V{q,t) = {(^T)r{q) - {(^H)r{q) = (-D- l)(cr|^)r(g) on r = ^. 

We use the following expansion for the mean curvature (compare with Eq. (4.12) of 


page 647 of mi and page 12 of [21], where instead of ( |4.3| ), the mean curvature of a 
{d- l)-sphere radius ro is defined to be ^): 


K{r = w) 


d-1 d-l 




q 


so that the linearisation of (4.14e) about r = q is 

13'j d-l 


{f^T)r{q) + W{q,t) = J^'j-XAi^T)r{q) +V{q,t)) 


on r = q. 


Finally, by the relation Vf (\x\) ■ i/ = f'(r) for a; e S, we have that 

V{VZ)-u\r=q = dr{V{r,t)Z{e,^))\r=g = Vr{q,t)Z{B,^), 

and so we obtain the following system for the perturbations U, V, W and 6 from linearising 
(4.14) about the unperturbed interface r = q: 

d \ d 

Wrr + = -VV 

Vrr+'^Vr+^V = CV 

Urr + '^Ur + = 0 

U-V = {D-l){aH)r{q) on r = q, 


{Vt + X^crT)r{q) + W + x^V ■■ 


/37 d - 1 




2 

( ct ^ - DaH)rr{q) = DUr - Vr 

= -(l^T)rrO 


We complete (4.18) with the following boundary conditions; 

VF(r = 0, t) < oo, E(r = 0 ,t)<oo, U{r = R,t) = 0- 


in r < q, 

(4.18a) 

in r < q, 

(4.18b) 

in r > q, 

(4.18c) 

on r = q, 

(4.18d) 

on r = q, 

(4.18e) 

on r = q, 

(4.18f) 

on r = q. 

(4.18g) 

= 0. 

(4.19) 
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4.3 Solutions to the perturbed system 

Recalling the definition of in (4.17), we see that the general solution for (4.18c) is 


TT( ^ forfi=2, 

’ I Fo(t)r^ + for d = 3. 


(4.20) 


We observe that the ODE (4.18b|) in d = 2 is a scaled modified Bessel’s equation (see 


(4.8)), while (4.18b) in d = 3 is a scaled modified spherical Bessel’s equation. Due to the 


boundary condition (4.19), we see that the general solution to (4.18b) is given by 


E(r,t) = 


{F 2 {t)Ii{Ar) ford = 2, 
|F 2 (t)i;(Ar) for d = 3, 


(4.21) 


where Ii{x) is the modified Bessel function of the first kind, dehned in (4.7), while ii{x) 
is the modified spherical Bessel function of first kind that satisfies 

- ^H{r) = ii{r), ^^(0) < oo VZ > 0. 

dr^ r dr r^ 


Again, due to the boundary condition (4.19 ) 2 , V(r, t) does not contain any terms involving 


the modified spherical Bessel function of the second kind. 


For (4.18a), we see that VE is a sum of the solution to the homog eneou s equation 
(4.18c) and the particular solution Due to the boundary condition (4.19) for W, we 


find that the general solution to (4.18a) is 


J 


Wir,t) = Fsity - - Vir,t). 


(4.22) 


With these solutions (4.20), (4.21), and (4.22), we use the relations (4.18d), (4.18e), (4.18f) 


in order to simplify the resulting differential equation (4.18g) for 6. Let 


Q(A,q) := 


A( 7 cosh(Ag) - sinh(Ag) 




then from (4.18d), (4.18e), and (4.18f) we obtain the following relations: 

for d= 2, 


Foq^ + - F2li{Aq) = {D-l)f 


Fog' + Fig-'-i - F2ii{Aq) = {1 - D)f 


(4.23a) 


for d = 3, 


(x<p - f) {b2AIi{Aq) + F2li{Aq)) + f g + Fag' = 


for d=2, 


{x^ - f) {b3Q{A,q) + F2ii{Aq)) + fq + F^ = for d = 3, 


(4.23b) 


'Cb 2 h{Aq) = F(/Fog'-i - Fi/g-'-^) - F 2 AI/(Ag) for d = 2, 

(4.23c) 

= Z4(/Fog'-i - F^ + l)g-'-2) - F 2 Az((Ag) for d = 3, 
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where we have used that 

+ Xv<yT)r{q) 


^0' +(xip - f )^^ 2 -fi(Ag) ford = 2, 
i(l + {Xip-^)h{t)Q{Kq) ford = 3, 


(4.24) 


and by (4.4d)i, 


(o-f - DaH)"{q) = Ca^iq) - ^^—^(cr^Yiq) + = Ca^iq)- 


(4.25) 


Also, from (4.4d)i, we observe that the following relations hold 

, , {b 2 Ah(Aq) = D^ ford = 2 , 

(cjy)'(g) = L»(a|^)'(g) ^ i i U y; 9 
KTJK^J ford = 3 . 

Together with the relation 

{fiY)"{q) + VaYiq) =A- ^^(/xy)'(g) 

q 


(4.26) 


= A- 




'i + ^Df ford =2, 
ford=3, 


and the relations (4.23), we can simplify (4.18g) in order to obtain the following differential 

1{1 + 2){1-1) 


equation for the perturbation size 6: 

-(/-l)-^(/x^-(U2d))|)-/37^ 


_ A, 

d dt 3 


2g3 


W"' [ix^ + KD - 1)|) + ^ _ (/ + j for d 


(4.27) 


and 


1 ^ A 


[ix^ + KD - 1 )|) + A _ (/ + iD)^^ for d = 2. 


2g3 


(4.28) 


Consequently, using (4.11), and (4.26), we obtain the following differential equations for 
the shape perturbation 


gd/d\_ldd Idg 
d dt \g/ d dt g dt 


A 02 

qz 


(^X¥>-(^ + 2 T>)^) -^7 


2g3 


(/X^ + - 1)|) + ) 


for d = 2, 


(4.29) 


/^ _ 03 

qo 


(^X<p-(^ + 3T>)^) -/37 


K/ + 2 )(Z- 1 ) 


+ A’od'-' [ix^ + KD-l)^y^ (/X^ - (/ + /D + dd)^) 


for d = 3. 
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Finally, we mention that the time-dependent constants Fq and -Fi can be computed as 
follows: Due to (4.19)3, we have 


ford = 2 , 
\-FiR-^^-^ ford = 3. 


Moreover, by (4.23a) and (4.23c), we obtain 


Cb2lo{Aq) + Af^{l-D)^ 

Il{Aq) q 

j, (Di^ m 

H Ii{Aq)W 

Cbs smhjAq) ^ ^^ i'ljAq) ^^03 

q ii{Aq)^ q^ 


= -Fi 


IDq^-^ {l + l)D A{Aq) 

d22^+i ql+2 



for d = 2 , 




for d = 3, 


(4.30) 


(4.31a) 


(4.31b) 


respectively. 

We observe that the active transport parameter A enters into the radial solutions (4.6), 
the differential equations (4.11), (4.27), (4.28), and (4.29) only via the time-dependent 
constants 02,03, 62 , ^3, C2 and C3. 


4.4 Effect of active transport on linear stability 

We now investigate the effect of active transport on the linear stability of the system. To 
compare with [TTj, we consider the choices 


C = 1, A = 1 , CTcx, = 1, 

and neglect Fq in ( 4.29[ ). This implies that (4.31) becomes 


C 63 sinh(Ag) 


q 


ii{Aq) 


q" 


-M 


q"- Il{Aq) q^ 

{Ul)D ^ i\{Aq) 1 




+ A 


il{Aq) 


for d = 2 , 


for d = 3. 


We define 
02 


qh{q) 


as = 


DIo{q) - qlog{qlR)Ii{q) ’ 
Rq{l - gcoth(g)) 


b2 
h = 


D 


DIo{q) - qlog{qlR)Ii{q) ’ 
DRq 


{R-q){qcotli{q) - 1) + DR ' (72 - g)(gcosh(g) - sinh(g)) + iAi?sinh(g) ’ 

so that 02 = 02(1 + 2A), 62 = 62(1 + 2A), 03 = 03(1 + 2A) and 63 = 63(1 + 2A), where 02 , 03 , 
62 and 63 are as defined in (4.10). A short computation yields that 

-Da^ _ A)(coth(g)-p 


02-.= 


Da2 


D/i(g)//o(g) 


q D-qlog{qlR)Ii{q)IIo{q) 


Cs:= 


g 2 ^ + ■ 


Using the following relations for the modified Bessel functions and modified spherical 
Bessel functions of the first kind: 

f I f I j TT 

Ii{z) =-Ii{z) + Ii+i{z), i[{z) = -ii{z) + ii+i{z), ii{z) = \ —Ii^iiz), 

Z Z y ZZ 2 
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and the relations 


hioiq) 

— sinh(g) 
63- 


DIo{q) 


Hq) 


DIi{q) 


DIo{q) -qlog{q/R)Ii{q) h{q) DIq - qlog{q/R)Ii{q) 

_ im __ C3 

{R-q){qcoih.{q)-1) + DR coth(q') - ^ ’ 


^^2 


Hq) 

h{qy 


we find that 


^1 = 




( h(q) . (I , li+i{q) 


D 


q ii(q) 


)) 


-{l + 2X)y^‘^C3 


( 1 , l:zll ( h+a/al'?) , 1 

\coth(q)-l/q D \Ii+i/2(q) q 


for d = 2, 


)) 


(4.32) 




for d = 3. 


Substituting Fq = 0 and A = 0 in (4.32) and (4.29), we obtain the differential equation for 
the shape perturbation as derived in Eq. (89) of [H] with the notation := \p^- 


d 5 


Next, we find, for given V, D, Xtp: arid /3, a critical value Ale such that = 0, i.e., 
the shape perturbation (^) is a constant. This critical value Ale is given by 


jAr 


Px^^-^ + {1 + 2X)2C2-^ ' 


Dq 


+ il + 2X)2C2ix^-il + D)V) 


( io(q) . ]-D (I . ii+i(q) \\ 
[h(q) + g U h(g) )) 


for d = 2, 


and 


Ale = PX 


3{l + 2){l-l) 
2q^ 


+ (1 + 2A)3C*3 


(1 + 


^-r)'P-x, 

Dq 


+ {l + 2X)3C3{x,2-{l + D-f)V) 


( 1 , UAl ( h+3/2('?) . l \\ 

\coth(q)-l/q D \Ii+i/2(q) qj) 



+ I + 


h+3/2(g) \ 
^^+ 1 / 2 ( 9 )) 


for d = 3. 


We point out that, when A = 0, the expression for Ale coincides with Eq. (90) of [TTj with 
Q~^ '■= G~^t = I/ 37 . We now look at Ale as a function of q for the following parameter 
values: 


= 0.05, ^ = 0.1, D = l, 1 = 2, R = 13. 

With these choices, we obtain 

+ (1 + 2A)(x^ - 0 . 2)202 U - M for d = 2, 

A = ^ 

“,3C3(l + 2A)(x,(i'-i) + (A-0.15>')) for<i = 3. 


(4.33) 
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where 


^ _ h{q)IIo{q) (coth(g)-^) 

^ l-glog(g/i?)/i(g)//o(g) ’ ^ i + g^(coth(g)-i) ’ 

Y - ^ -^o(g) y _ 1_1 

A + q^-^h{q)' 5 + Q ^V5(g) coth(g) -1/g ' 

^ ^ 1512(g) 

Numerically, we find that C 2 , C 3 , X and Y are positive for q e (0,13]. Moreover, 

X--<0, y--<0, —-0.15y>0 Vge(0,13]. (4.34) 

q q 4:q 

We note that A is the apoptosis parameter and Ac divides the phase portrait into regions 


A as a ftinrtw® of q {br rf = 2. \. = 0 and A = {0. 2.5.10} A as a fiuirtkin of 9 for d = 3. = 0 ajid A = {0. 2. r>. 10} 




(a) 


(b) 


/t as a AmrtKHi of <7 ftir rf = 2. = {0.1,0.3} and A = {2.4} /I as a fiuirtimi of g Inr d = 3. = {0.3.1.7} and A = {2.4} 




(c) (d) 


Figure 1: Effects of A on the critical apoptosis parameter Ac as a function of the unper¬ 
turbed radius q in 2d and 3d with /?7 = 0.1, V = 0.1, Z1 = 1, f = 2, i? = 13. 


of stable growth for low apoptosis (the region A < Ac) and regions of unstable growth for 
high apoptosis (the region A > Ac) for a given mode 1. Thus, from (4.33) and (4.34), we 
observed the following: 

1. In the absence of chemotaxis, = 0, increasing A will increase the value of Ac- From 
Figures |l(a) and |l(b)[ the curves are pushed upwards, and so the region of stable 
growth for low apoptosis is enlarged. In particular, active transport has a stabilising 
effect on the perturbations in the absence of chemotaxis. 
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2. In dimension d = 2, while x<p < 0-2, active transport has a stabilising effect on the 
perturbations. When Xip > 0.2, the perturbations are now amplified by the presence 
of active transport. In Figure 1(c), we see that, as A increases, the curves are pushed 
up for x<fi = 0.1, while the curves are pulled down for x<p = 0.3. Similarly, in dimension 
d = 3, we find that 

0.25/,-(USn,) ^ ^ ^ ^ ^ 

and from Figure l(d)[ we see that, as A increases, the curves are pushed up for 
= 0.3, while the curves are pulled down for = 1.7. 


5 Numerical Computations 


In this section we first derive a finite element approximation of (3.63) and then we display 
some numerical results obtained using this approximation. We concentrate on (3.63), 
however approximations of other variations of the model follow in a natural way. In the 
approximation we take '!'(¥?) to be the double obstacle potential given in (|3.1[ ). This choice 
of potential leads to (3.63b) taking the form of a variational inequality (|3.4|). 


Finite element approximation 


Let T be a regular triangulation of Q into disjoint open simplices, associated with T is 
the piecewise linear finite element space 


Sh ■■= {if e C'°(L!)| e Pi(r) V T e r} c 


where we denote by Pi{T) the set of all affine linear functions on T. We now introduce 
a finite element approximation of (3.63) in which we have taken homogeneous Neumann 
boundary conditions for p and /i, and the Dirichlet boundary condition cr = cjs € M on 512: 
Find 


p^^h^K^,:={x€Sh\\x\<l}, a'fl€S^-.= {x€Sh\x = ^Bondn} 

such that for all r]h e Sh, (h e Kh and Xh e 5° := {x^ Sh\ X = 0 on 512}, 

-{‘Ph - ‘Ph~^,r]h)h + = {iVah~^ - A){<pI + l),Vh)h, (5.1a) 

r 

(5.1b) 

\ S /h 

{V{ipl)Val, VXh)h - VXh)h = + l),Xh)h, (5.1c) 

where m{(p) = ^(1 + y?)^, r denotes the time step, (r/i,?? 2 ) denotes the L? inner product 
and {r]i,r] 2 )h •= dx where on each triangle TTh is taken to be an affine 

interpolation of the values of r]ir ]2 at the nodes of the triangle. 

We note that since the interfacial thickness is proportional to e in order to resolve the 
interfacial layer we need to choose h « e, see [18] for details. Away from the interface h 
can be chosen larger and hence adaptivity in space can heavily speed up computations. 
In fact we use the finite element toolbox Alberta 2.0, see |l2j, for adaptivity and we 
implemented the same mesh refinement strategy as in [5|, i.e., a fine mesh is constructed 
where < 1 with a coarser mesh present in the bulk regions = 1 - 

We begin our numerical results by following the authors in m in comparing solutions 
obtained from a simplified form of the diffuse interface model with exact solutions to a 
sharp interface limit model. 
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5.1 Comparison with a sharp interface limit solution 

In Figures and we display results obtained from the growing circle tumour test case 
introduced in Section 4.2 of [30]. To this end we consider the simplified model on a circular 
domain Q with radius R: 


14^/2 2^ 

dtp = A/x + P )<7, 

(5.2a) 

e TT 

p = -^'{p) -eAp, 
e 

(5.2b) 

O-Ag ^^^(1 p'^)g. 

(5.2c) 


e TT 


Here (p and p satisfy homogeneous Neumann boundary conditions, and a satishes the 
Dirichlet boundary condition a = gr e M on 912. We take the radially symmetric case of 
an initial circular tumour with initial radius 0.25. From m we have that the solution to 


the sharp interface limit of (5.2) is given by 


a{r,t) = 


^p{t) 




r < p{t), 


(5.3) 


where ^ 2s/2p{t)\og{p{t) jR) ’ ^ being constant and p{t), which is the radius of the 

tumour, is determined by numerically solving the ODE p'{t) = \/ 2 a{p{t),t) with initial 
condition p(0) = 0.25. 

We set 12 = 10 and gr = 2, however for the diffuse interface computations we did not 
solve the problem in the whole of 12 instead we solved it on a circular domain with radius 2 
with the time dependent Dirichlet boundary condition G{x,t) - o-£)(|x| ,t) computed from 
) with r = 1. We set r = 1.0e“^, the minimal diameter of an element hmin = 7.8125-10“^ 
and the maximal diameter hmax = 3.125 • 10“^. 

In Figurej^we display the diffuse interface solutions p and cr at t = 0, 0.2, 0.4 obtained 
with e = 0.05. In the plots of (p we include the sharp interface limit solution of the tumour 
position. In Figure]^ we examine the convergence of the diffuse interface solution to the 
sharp interface limit solution as e tends to zero. In Figure |^a) we plot the radius of the 
growing tumour for the diffuse interface model with e = 0.1, 0.075, 0.05 together with the 
sharp interface limit solution p{t). In Figure j^b) we plot the solution g of the diffuse 
interface model with e = 0.1, 0.075, 0.05 together with the sharp interface limit solution 
G at t = 0.1. From this figure we see that as e decreases the diffuse interface solution 
converges to the sharp interface limit solution. 



5.2 Solutions of (5.1) 


We now investigate the influence of the parameters V, Xtp S'lid A in Model (3.63). In all 
computations we set f2 = (-12.5,12.5)^, M = 0, D = 1, /3 = 0.1, C = 2, ub = 1, r = 1.0e“^, 
the minimal diameter of an element hmin = 4.888 ■ 10“"^ and the maximal diameter hmax = 
5 ■ 10“^. Unless otherwise specified we take e = 0.01. 


Influence of the proliferation rate V 

In Figures]^ andwe investigate the influence of V. We set x^p = ^ a-ncl A = 0.03. In Figure 
[^we set P = 0.5 while in Figure]^ we set V = 0.1, and in both sets of figures we display ip 
(top row) and g (bottom row) at times t = 5, 10,13. From this figure we see taking V = 0.5 
gives rise to Angers that are thicker than the ones resulting from V = 0.1. 
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Figure 2: Approximate solutions of (5.2) at t = 0 (left), t = 0.2 (centre) and t = 0.4, ip (top 
row), a bottom row. The black line in the ip solutions denotes the corresponding sharp 
interface solution. 




(a) radius versus time 


(b) a at t = 0.1 


Figure 3: Comparison of diffuse interface model (5.2) with the sharp interface solution. 


Influence of the chemotaxis parameter Xip 

In Figures]^ andwe investigate the influence of X(fi- We set V = 0.1 and A = 0. In Figure 
[^we set Xifi = ^ while in Figure]^ we set Xif = 10, and in both sets of figures we display ip 
(top row) and a (bottom row). The results for Xy, = 5 are displayed at times t = 5, 10,20, 
while the results for Xtp = 10 are displayed at times t = 2.5,5,10. From these figures we 
see that, akin to the results in [H], for both values of Xip after some time fingers develop, 
and thereby increasing the surface area of the tumour to allow for better access to the 
nutrient. For the larger value of Xv th® formation and evolution of the fingers is quicker 
and the fingers are slimmer. 
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0.0607 


0.0295 


0.000342 





Figure 4: Solutions of (5.1) with A = 0.03, Xv? = 5, P = 0.5, at t = 5,10,13. 


Influence of the active transport parameter A 


In Figures]^- 10 we investigate the influence of A. We set V 
we show if (top row) and a (bottom row) at t 


0.1 and x^p 
0 (left), A = 


5. In Figure 
0.07 (centre) 


4, with A 

and A = 0.09 (right). From this figure we see that when A = 0 the variation of a across 
the interfacial region is smooth while taking A > 0 leads to a drastic change in a. This 
change in a can be seen better in Figure where we show plots of and a along a line 
that spans the interfacial region. The scales for a and ip are shown on the left and right 
axes respectively. 

Here we see that the change in a across the interfacial region is more pronounced for 
larger values of A. In Figure [T0| we display the influence of e on the change in a across the 
interfacial region, we set A = 0.07 and plot a along a line that spans the interfacial region 
for £ = 0.04, 0.02, 0.01. From this figure we see the convergence of u as e decreases. In 
Figure 10 the jump in a across the interfacial region for e = 0.01 is 0.1327 ~ 2A which is 


consistent with the formal asymptotic analysis, recall (3.65). 


5.3 Numerical computations with Darcy flow 

For positive constants mo and K, we now consider the model 


divr? = oF, 

(5.4a) 

V = -K{Vp-{p + XpCF)Vip), 

(5.4b) 

dfip + div (ipv) = moA|U + p^F, 

(5.4c) 

p = -^'{ip) - (3£Aip - Xpcr, 

£ 

(5.4d) 

0 = div {'D{(p){Va - AV(^)) - -Ca{ip + 1), 

(5.4e) 
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Figure 5: Solutions of (5.1) with A = 0.03, Xi/? = 5, P = 0.1 at t = 5,10,13. 


where we recall that a -==— nc ;= F 

(3.62). As additional boundary condition we prescribe 

p = 0 on 


\{V(J - + 1), and V is defined in 


while we take homogeneous Neumann boundary conditions for (p and //, and the Dirichlet 
boundary condition fi = (Tb e M on dO.. Recalling the finite element spaces Kh, Sh, Sj^ and 
defined at the start of Section for the double-obstacle potential 
the following scheme for the above system: Find 

such that for all Sh, Ch ^ Kh and x/j e 5°, 

- ^h~^,Vh)h + rnoiVpl, Vvh) 

T 

= y((p<"‘ -.4)w++ !),%)„ 


+ K{Vpt^ ■ Vpt' - {pt' + xX“') \^‘fVf,Vh)h, (5.5a) 

if^h + ~^h^'^X<p<^h^,Ch-^h') ^V(C/i - (5.5b) 

V ^ / h 

{V{pl)Val, VXh)h - X{V{pl)Vpl, VXh)h = -]f{al{pl + l),Xh)h, (5.5c) 

i^Ph^X/Xh) = HpI + Xv^D^^h^ X/Xh)h + i^ii'Pcrh - + l),Xh)h- (5.5d) 


As initial condition for p and p, we always choose = 0 and = 0. We perform three 
different numerical simulations in which we vary the tumour and healthy cell densities. 
The three cases are given as follows: 


(3.1), we propose 
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Figure 6: Solutions of (5.1) with V = 0.1, A = 0, Xyj = 5 at t = 5,10,20. 


(Case (1)) a = 0 and ps = with = p 2 = 1 so that we solve for 

divn = 0, dtp + S/p ■ v = moAp + {Va - A){p + 1); 


(Case (2)) a = I and ps = ‘^ with pi = p 2 = f so that we solve for 

divr; = -{Va - A){p + 1), dtp + div (pv) = moAp + [Va - Al)((/? + 1); 
3 

(Case (3)) a = -| and ps = ‘^ with Pi = |, /O 2 = i so that we solve for 


divn = — {Va - Al)(</? + 1), dtp + div {pv) = moAp + {Va - A){p + 1). 

3 

We always take cts = 1, /3 = 0.1, V = 0.1, Al = 0, C = 1, = 10, e = 0.01, mo = 1 K = 0.01, 

A = 0.03 and D = 1. 


In Figure 11 we display the solutions of the Darcy flow model (5.5) for case (1) at 
t = 1.5; the left plot is of p, the centre plot of a and the right plot is of p. In the left 


plot of Figure 12 we display a zoomed in plot of p at t = 1.5 obtained from the Darcy 


flow model (5.5) for case (2) together with the zero level line of p't^{x) (depicted in black) 
with K = a = 0 (which is equivalent to (5.1)). In the right plot we show the influence of 
a on the position of the tumour, we display a zoomed in plot of the solution y? at t = 1.5; 
the black, white, and blue lines are the zero level lines of p for case (1), case (2) and 
case (3), respectively. One observes that the model variant with Darcy flow enhances the 
growth velocities of the tumour. In addition, the velocity is largest when the density of 
the tumour is smaller than the density of the healthy cells. 
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Figure 7: Solutions of (5.1) with V = 0.1, A = 0, Xv? = 10 at t = 2.5,5,10. 
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Figure 9: Comparison of a across the interfacial region for (5.1) at t 
A = 0.07 (centre) and A = 0.09 (right). 


4 with A = 0 (left), 



Figure 10: Convergence of a as e decreases for (5.1) with A = 0.07 at t = 4. 
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line denoting the zero level line of (5.1). The right plot displays the zero level 

lines of pi^{x) from (5.5) for case (1) (black line), case (2) (white line) and case (3) (blue 
line) at t = 1.5. 
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